In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import Image
import plotly.express as px
import os
import plotly.graph_objects as go
In [2]:
act = pd.read_csv('ActualEvapoTranspiration_PMFAO_4.csv',skiprows=6, sep=',', parse_dates=[0], na_values=-9999,usecols=[1,2])
act.columns = ['Datetime','EvapoTranspiration']
pot = pd.read_csv('PotentialEvapoTranspiration_PMFAO_4.csv',skiprows=6, sep=',', parse_dates=[0], na_values=-9999,usecols=[1,2])
pot.columns = ['Datetime','EvapoTranspiration']
wc = pd.read_csv('wcEvapoTranspiration_PMFAO_4.csv',skiprows=6, sep=',', parse_dates=[0], na_values=-9999,usecols=[1,2])
wc.columns = ['Datetime','EvapoTranspiration']

fig = px.line()
fig.add_trace(go.Scatter(x=pot['Datetime'], y=pot['EvapoTranspiration'], mode='lines', name='Potential ET'))
fig.add_trace(go.Scatter(x=act['Datetime'], y=act['EvapoTranspiration'], mode='lines', name='Actual ET'))
fig.add_trace(go.Scatter(x=wc['Datetime'], y=wc['EvapoTranspiration'], mode='lines', name='Water Content ET'))

fig.update_layout(
    title='Penman-Monteith EvapoTranspiration with different soil moisture values',
    xaxis_title="Time",
    yaxis_title="EvapoTranspiration [mm]",)
fig.update_layout(yaxis_range=[-0.1, 1.1])
fig.update_xaxes(rangeslider_visible=True)
fig.show()
In [3]:
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc
import matplotlib.style as style 
import matplotlib.colors as cl
import matplotlib.dates as mdates
import matplotlib
from geoframepy.timeseries import io_csv
%matplotlib inline
my_cmap = matplotlib.colors.LinearSegmentedColormap.from_list("", ["#ece7f2","#d0d1e6","#a6bddb","#74a9cf","#3690c0","#0570b0","#045a8d","#023858"])
colors = ['#a6cee3', '#1f78b4', '#b2df8a', '#33a02c', '#fb9a99', '#e31a1c', '#fdbf6f', '#ff7f00']
folderPath = "C:/Users/E580/Desktop/UNI/Idrologia/Esercitazione3"
In [4]:
rain = io_csv.pandas_read_OMS_timeseries(folderPath + '\Precout_long.csv',datetime_index=False)
rain.head()
Out[4]:
Datetime 0
0 2021-07-05 00:00:00 0.0
1 2021-07-05 01:00:00 0.0
2 2021-07-05 02:00:00 0.0
3 2021-07-05 03:00:00 0.0
4 2021-07-05 04:00:00 0.0
In [5]:
ds = xr.open_dataset(folderPath + '\pistachio2_output_long_0000.nc')
ds
Out[5]:
<xarray.Dataset>
Dimensions:           (depth: 600, dualDepth: 600, time: 481)
Coordinates:
  * depth             (depth) float64 -2.998 -2.993 -2.988 ... -0.0075 -0.0025
  * dualDepth         (dualDepth) float64 -3.0 -2.995 -2.99 ... -0.01 -0.005
  * time              (time) datetime64[ns] 2021-07-05 ... 2021-07-25
Data variables:
    psi               (time, depth) float64 ...
    psiIC             (depth) float64 -0.0025 -0.0075 -0.0125 ... -2.992 -2.998
    T                 (depth) float64 273.1 273.1 273.1 ... 273.1 273.1 273.1
    theta             (time, depth) float64 ...
    saturationDegree  (time, depth) float64 ...
    darcyVelocity     (time, dualDepth) float64 ...
    waterVolume       (time, depth) float64 ...
    error             (time) float64 -1.926e-15 -2.006e-15 ... -2.107e-15
    topBC             (time) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    bottomBC          (time) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    runOff            (time) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    controlVolume     (depth) float64 0.005 0.005 0.005 ... 0.005 0.005 0.005
Attributes:
    Description_of_the_problem:         WHETGEO-1D ex03 wet initial condition.
    Top_boundary_condition:             Top Coupled
    Bottom_boundary_condition:          Bottom Dirichlet
    path_top_boundary_condition:        C:/Users/E580/Desktop/UNI/Idrologia/E...
    path_bottom_boundary_condition:     C:/Users/E580/Desktop/UNI/Idrologia/E...
    path_grid:                          C:/Users/E580/Desktop/UNI/Idrologia/E...
    time_delta:                          
    swrc_model:                         {Water Depth,Van Genuchten}
    soil_hydraulic_conductivity_model:  {null,Mualem Van Genuchten}
    interface_conductivity_model:       max
xarray.Dataset
    • depth: 600
    • dualDepth: 600
    • time: 481
    • depth
      (depth)
      float64
      -2.998 -2.993 ... -0.0075 -0.0025
      units :
      m
      long_name :
      Soil depth.
      array([-2.9975e+00, -2.9925e+00, -2.9875e+00, ..., -1.2500e-02, -7.5000e-03,
       -2.5000e-03])
    • dualDepth
      (dualDepth)
      float64
      -3.0 -2.995 -2.99 ... -0.01 -0.005
      units :
      m
      long_name :
      Dual soil depth.
      array([-3.   , -2.995, -2.99 , ..., -0.015, -0.01 , -0.005])
    • time
      (time)
      datetime64[ns]
      2021-07-05 ... 2021-07-25
      long_name :
      Time.
      array(['2021-07-05T00:00:00.000000000', '2021-07-05T01:00:00.000000000',
       '2021-07-05T02:00:00.000000000', ..., '2021-07-24T22:00:00.000000000',
       '2021-07-24T23:00:00.000000000', '2021-07-25T00:00:00.000000000'],
      dtype='datetime64[ns]')
    • psi
      (time, depth)
      float64
      ...
      units :
      m
      long_name :
      Water suction.
      [288600 values with dtype=float64]
    • psiIC
      (depth)
      float64
      ...
      units :
      m
      long_name :
      Initial condition for water suction.
      array([-2.5000e-03, -7.5000e-03, -1.2500e-02, ..., -2.9875e+00, -2.9925e+00,
       -2.9975e+00])
    • T
      (depth)
      float64
      ...
      units :
      K
      long_name :
      Temperature.
      array([273.15, 273.15, 273.15, ..., 273.15, 273.15, 273.15])
    • theta
      (time, depth)
      float64
      ...
      units :
      long_name :
      theta for within soil and water depth.
      [288600 values with dtype=float64]
    • saturationDegree
      (time, depth)
      float64
      ...
      units :
      1
      long_name :
      saturation degree.
      [288600 values with dtype=float64]
    • darcyVelocity
      (time, dualDepth)
      float64
      ...
      units :
      m/s
      long_name :
      Darcy velocity.
      [288600 values with dtype=float64]
    • waterVolume
      (time, depth)
      float64
      ...
      units :
      m
      long_name :
      Water volume in each control volume
      [288600 values with dtype=float64]
    • error
      (time)
      float64
      ...
      units :
      m
      long_name :
      Volume error at each time step.
      array([-1.925868e-15, -2.006345e-15, -1.901498e-15, ..., -1.890581e-15,
       -2.222494e-15, -2.107430e-15])
    • topBC
      (time)
      float64
      ...
      units :
      mm
      long_name :
      Rainfall heights
      array([0., 0., 0., ..., 0., 0., 0.])
    • bottomBC
      (time)
      float64
      ...
      units :
      long_name :
      array([0., 0., 0., ..., 0., 0., 0.])
    • runOff
      (time)
      float64
      ...
      units :
      m/s
      long_name :
      run off
      array([0., 0., 0., ..., 0., 0., 0.])
    • controlVolume
      (depth)
      float64
      ...
      units :
      m
      long_name :
      dimension of each control volumes
      array([0.005, 0.005, 0.005, ..., 0.005, 0.005, 0.005])
  • Description_of_the_problem :
    WHETGEO-1D ex03 wet initial condition.
    Top_boundary_condition :
    Top Coupled
    Bottom_boundary_condition :
    Bottom Dirichlet
    path_top_boundary_condition :
    C:/Users/E580/Desktop/UNI/Idrologia/Esercitazione3/Precout_long.csv
    path_bottom_boundary_condition :
    C:/Users/E580/Desktop/UNI/Idrologia/Esercitazione3/ex03_dirichlet_wet_long.csv
    path_grid :
    C:/Users/E580/Desktop/UNI/Idrologia/Esercitazione3/pistachio2.nc
    time_delta :
    swrc_model :
    {Water Depth,Van Genuchten}
    soil_hydraulic_conductivity_model :
    {null,Mualem Van Genuchten}
    interface_conductivity_model :
    max

plot dei grafici¶

In [19]:
fig, ax = plt.subplots(2, 1, sharex='col', figsize=(16, 10)) 
norm = matplotlib.colors.Normalize(vmin=np.amin(ds.psi.values), vmax=1.0)

ax[0].plot(ds.time.values, rain['0'])
ax[0].set_ylabel('[mm/s]')

ax1=ax[0].twinx()
ds.waterVolume.sel(depth=0.0, method='nearest').plot(ax=ax1,color='0.3')
ax1.set_ylabel('Water ponding [mm]')
ax1.set_title(None)


a = ds.psi.plot(ax=ax[1], y='depth',add_colorbar=False, cmap=my_cmap, norm=norm)#, cbar_kwargs={'label': '$\psi$ [m]'})

ax[1].set_ylabel('Depth [m]')
ax[1].set_xlabel('Time [dd hh:mm]')


ax[1].xaxis.set_major_formatter(mdates.DateFormatter('%d %H:%M'))
ax[1].xaxis.set_tick_params(rotation=45)

cbar_ax = fig.add_axes([0.95, 0.19, 0.02, 0.36])
sm = plt.cm.ScalarMappable(cmap=my_cmap, norm=norm)
cbar = fig.colorbar(sm, cax=cbar_ax)
cbar.set_label('$\\psi$ [m]', rotation=0, labelpad=22)

plt.show()
In [20]:
fig, (ax1, ax2) = plt.subplots(1, 2, sharey='row', figsize=(16, 10)) 

ds.psiIC.plot(ax=ax1, y='depth', label='Initial condition', color='0.3')
ds.psi.sel(time='2021-07-05 00:00').plot(ax=ax1, y='depth', label='2021-07-05 00:00', color=colors[0])
ds.psi.sel(time='2021-07-08 00:00').plot(ax=ax1, y='depth', label='2021-07-08 00:00', color=colors[1])
ds.psi.sel(time='2021-07-10 00:00').plot(ax=ax1, y='depth', label='2021-07-10 00:00', color=colors[2])
ds.psi.sel(time='2021-07-15 00:00').plot(ax=ax1, y='depth', label='2021-07-15 00:00', color=colors[3])
ds.psi.sel(time='2021-07-20 00:00').plot(ax=ax1, y='depth', label='2021-07-20 00:00', color=colors[4])
ds.psi.sel(time='2021-07-25 00:00').plot(ax=ax1, y='depth', label='2021-07-25 00:00', color=colors[5])

ax1.set_title('Water pressure \nat different times')

ax1.set_ylabel('Depht [m]')
ax1.set_xlabel('$\psi$ [m]')

plt.plot(np.nan, np.nan, label='Initial condition', color='0.3')
ds.theta.sel(time='2021-07-05 00:00').plot(ax=ax2, y='depth', label='2021-07-05 00:00', color=colors[0])
ds.theta.sel(time='2021-07-08 00:00').plot(ax=ax2, y='depth', label='2021-07-08 00:00', color=colors[1])
ds.theta.sel(time='2021-07-10 00:00').plot(ax=ax2, y='depth', label='2021-07-10 00:00', color=colors[2])
ds.theta.sel(time='2021-07-15 00:00').plot(ax=ax2, y='depth', label='2021-07-15 00:00', color=colors[3])
ds.theta.sel(time='2021-07-20 00:00').plot(ax=ax2, y='depth', label='2021-07-20 00:00', color=colors[4])
ds.theta.sel(time='2021-07-25 00:00').plot(ax=ax2, y='depth', label='2021-07-25 00:00', color=colors[5])

ax2.set_title('Water content \nat different times')

ax2.set_ylabel('')
ax2.set_xlabel('$\\theta$ [-]')


ax2.legend(bbox_to_anchor=(1.01,0.65), ncol=1)

plt.tight_layout()
plt.show()
In [21]:
act = pd.read_csv('LongActualEvapoTranspiration_PMFAO_4.csv',skiprows=6, sep=',', parse_dates=[0], na_values=-9999,usecols=[1,2])
act.columns = ['Datetime','EvapoTranspiration']
pot = pd.read_csv('LongPotentialEvapoTranspiration_PMFAO_4.csv',skiprows=6, sep=',', parse_dates=[0], na_values=-9999,usecols=[1,2])
pot.columns = ['Datetime','EvapoTranspiration']
wc = pd.read_csv('LongwcEvapoTranspiration_PMFAO_4.csv',skiprows=6, sep=',', parse_dates=[0], na_values=-9999,usecols=[1,2])
wc.columns = ['Datetime','EvapoTranspiration']

fig = px.line()
fig.add_trace(go.Scatter(x=pot['Datetime'], y=pot['EvapoTranspiration'], mode='lines', name='Potential ET'))
fig.add_trace(go.Scatter(x=act['Datetime'], y=act['EvapoTranspiration'], mode='lines', name='Actual ET'))
fig.add_trace(go.Scatter(x=wc['Datetime'], y=wc['EvapoTranspiration'], mode='lines', name='Water Content ET'))

fig.update_layout(
    title='Penman-Monteith EvapoTranspiration with different soil moisture values',
    xaxis_title="Time",
    yaxis_title="EvapoTranspiration [mm]",)
fig.update_layout(yaxis_range=[-0.1, 1.1])
fig.update_xaxes(rangeslider_visible=True)
fig.show()
In [ ]: